##########################################################################################

options(stringsAsFactors=F)
library(optparse)
library(patchwork)
library(ggplot2)
library(data.table)
library(tidyverse)
library(ggrepel)
library(grid)
library(maftools)
library(ggrepel)

##########################################################################################

option_list <- list(
  make_option(c("--gisticAllLesionsFile_IGC"), type = "character"),
  make_option(c("--gisticAmpGenesFile_IGC"), type = "character"),
  make_option(c("--gisticDelGenesFile_IGC"), type = "character"),
  make_option(c("--gisticScoresFile_IGC"), type = "character"),
  make_option(c("--all_data_by_genes_IGC"),type = "character"),
  make_option(c("--gisticAllLesionsFile_DGC"), type = "character"),
  make_option(c("--gisticAmpGenesFile_DGC"), type = "character"),
  make_option(c("--gisticDelGenesFile_DGC"), type = "character"),
  make_option(c("--gisticScoresFile_DGC"), type = "character"),
  make_option(c("--all_data_by_genes_DGC"),type = "character"),
  make_option(c("--dat_info_multi"),type = "character"),
  make_option(c("--dat_info_TCGA"),type = "character"),
  make_option(c("--out_path"), type = "character"), 
  make_option(c("--Class"), type = "character"),
  make_option(c("--geneprotein_list"), type = "character")
)

if(1!=1){
  
  gisticAllLesionsFile_IGC = 'all_lesions.conf_99.txt'
  gisticAmpGenesFile_IGC = 'amp_genes.conf_99.txt'
  gisticDelGenesFile_IGC = 'del_genes.conf_99.txt'
  gisticScoresFile_IGC = 'scores.gistic_IGC'
  
}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gisticAllLesionsFile_IGC <- opt$gisticAllLesionsFile_IGC
gisticAmpGenesFile_IGC <- opt$gisticAmpGenesFile_IGC
gisticDelGenesFile_IGC <- opt$gisticDelGenesFile_IGC
gisticScoresFile_IGC <- opt$gisticScoresFile_IGC
gisticAllLesionsFile_DGC <- opt$gisticAllLesionsFile_DGC
gisticAmpGenesFile_DGC <- opt$gisticAmpGenesFile_DGC
gisticDelGenesFile_DGC <- opt$gisticDelGenesFile_DGC
gisticScoresFile_DGC <- opt$gisticScoresFile_DGC
all_data_by_genes_IGC <- opt$all_data_by_genes_IGC
all_data_by_genes_DGC <- opt$all_data_by_genes_DGC
dat_info_multi <- opt$dat_info_multi
dat_info_TCGA <- opt$dat_info_TCGA
out_path <- opt$out_path
Class <- opt$Class
geneprotein_list <- opt$geneprotein_list

###########################################################################################
##设置的函数

gisticChromGet = function(gistic = gistic, fdrCutOff = fdrCutOff, ref.build = ref.build,class=class,sub=sub) {
  
  g = maftools::getCytobandSummary(gistic)
  g = g[qvalues < fdrCutOff]
  g[,Chromosome := sapply(strsplit(x = g$Wide_Peak_Limits, split = ':'), '[', 1)]
  g[,loc := sapply(strsplit(x = g$Wide_Peak_Limits, split = ':'), '[', 2)]
  g[,Start_Position := sapply(strsplit(x = g$loc, split = '-'), '[', 1)]
  g[,End_Position := sapply(strsplit(x = g$loc, split = '-'), '[', 2)]
  g.lin = maftools:::transformSegments(segmentedData = g[,.(Chromosome, Start_Position, End_Position, qvalues, Cytoband, Variant_Classification,Unique_Name)])
  
  g.lin <- subset(g.lin,Variant_Classification==sub)
  mb <- g.lin
  #mb$StartPos <- mb$Start_Position_updated
  #mb$EndPos <- mb$End_Position_updated
  
  gis.scores = maftools:::transformSegments(segmentedData = gistic@gis.scores, build = ref.build)
  gis.scores <- subset(gis.scores,Variant_Classification==sub)
  #gis.scores$amp = ifelse(test = gis.scores$Variant_Classification == 'Del', yes = -gis.scores$fdr, no = gis.scores$fdr)
  gis.scores$amp = gis.scores$G_Score
  
  fdrCutOff = -log10(fdrCutOff)
  
  ##设置主键
  data.table::setkey(x = gis.scores, Chromosome, Start_Position_updated, End_Position_updated)
  cyto_peaks_scores = data.table::foverlaps(y = gis.scores[,.(Chromosome, Start_Position_updated, End_Position_updated, amp,G_Score,fdr)],
                                            x = mb[,.(Chromosome, Start_Position_updated, End_Position_updated, Cytoband,qvalues,Unique_Name)])
  cyto_peaks_scores = cyto_peaks_scores[order(Cytoband, abs(amp), decreasing = TRUE)][Cytoband %in% mb$Cytoband][!duplicated(Cytoband)]
  cyto_peaks_scores = cyto_peaks_scores[complete.cases(cyto_peaks_scores)]
  cyto_peaks_scores$class <- class
  
  return(cyto_peaks_scores)
  
}


plotCytoband <- function(dat=dat,sig=sig,type=type){
p <- ggplot(data = dat, aes(x = Gscore_diff, y = logqvalues_IGC, color=note))+
  geom_point(size=0.5) +
  theme_classic()  +
  scale_color_manual(values = rev(c("#CC0000", "#000000"))) +
  xlab('GISTIC score difference') + 
  ylab(expression(paste(-log[10]," (FDR) IGC"))) + 
  geom_hline(yintercept=-log10(q_t),color='grey50',linetype="dashed") +
  geom_vline(xintercept=diff_t,color='grey50',linetype="dashed") +
  theme(legend.position="none",
        axis.text.x = element_text(color='black'),
        axis.text.y = element_text(color='black'),
        axis.ticks = element_line(color='black')) +
  #xlim(-1.5,1.5) + ylim(0,28) +
  # 添加标签：
  geom_text_repel(data = sig,
                  #min.segment.length = Inf,
                  max.overlaps = getOption("ggrepel.max.overlaps", default = 100),
                  aes(label = text),
                  color = '#CC0000',
                  point.padding = 0.3,
                  face = "italic" ,
                  #parse = TRUE,
                  size = 3,
                  hjust = 1,
                  vjust = 1
                  ) +
  #theme(text = element_text(face = "italic")) +
  labs(title = paste0("IGC vs DGC ","(",type,")"),size=2)

return(p)

}



trans <- function(num){
  up <- floor(log10(num))
  down <- round(num*10^(-up),2)
  text <- paste("P == ",down," %*% 10","^",up)
  return(text)
}



datWilcoxTest <- function(dat_gene1=dat_gene1,dat_gene2=dat_gene2,show_symbol=show_symbol){
  if(length(show_symbol) > 0){
  	##DUX2有重名的，all.gene.txt文件里面有两个DUX2|chr4和DUX2|chr10，需要指定一下
    for(i in 1:length(show_symbol)){
    	if(show_symbol[i]=="DUX2"){
    		show_symbol[i] <- "DUX2|chr4"
    	}else if(show_symbol[i]=="DUX4"){
    		show_symbol[i] <- "DUX4|chr4"
    	}else if(show_symbol[i]=="DUX4L6"){
    		show_symbol[i] <- "DUX4L6|chr4"
    	}else if(show_symbol[i]=="DUX4L5"){
    		show_symbol[i] <- "DUX4L5|chr4"
    	}else if(show_symbol[i]=="DUX4L3"){
    		show_symbol[i] <- "DUX4L3|chr4"
    	}else if(show_symbol[i]=="DUX4L2"){
    		show_symbol[i] <- "DUX4L2|chr4"
    	}
  	}
  dat_gene1 <- subset( dat_gene1 , Gene_Symbol %in% show_symbol )
  dat_gene2 <- subset( dat_gene2 , Gene_Symbol %in% show_symbol )
  
  rownames(dat_gene1) <- dat_gene1$Gene_Symbol
  rownames(dat_gene2) <- dat_gene2$Gene_Symbol
  result <- c()
  for(gene in show_symbol){
    print(gene)
    tmp1 <- subset(dat_gene1,Gene_Symbol==gene)
    tmp2 <- subset(dat_gene2,Gene_Symbol==gene)
    

    tmp1 <- data.frame( use_sample = colnames(tmp1[,-c(1:3)]) , copy_number = as.numeric(tmp1[,-c(1:3)]) , type = "IGC" )
    tmp2 <- data.frame( use_sample = colnames(tmp2[,-c(1:3)]) , copy_number = as.numeric(tmp2[,-c(1:3)]) , type = "DGC" )
    
    ## 同一患者存在多个的，取中位数
    tmp1 <- merge( tmp1 , dat_info[,c("Normal" , "use_sample")] , by = "use_sample" )
    tmp2 <- merge( tmp2 , dat_info[,c("Normal" , "use_sample")] , by = "use_sample" )
    tmp1 <- tmp1 %>%
      group_by( Normal , type  ) %>%
      summarize( copy_number = median(copy_number) )
    tmp2 <- tmp2 %>%
      group_by( Normal , type  ) %>%
      summarize( copy_number = median(copy_number) )    
    p <- wilcox.test( tmp1$copy_number , tmp2$copy_number )$p.value
    
    if( p < 0.01 ){
      p_text <- trans(p)
    }else{
      p_text <- paste0( "P == " , round(as.numeric(p) , 2) ) 
    }
    
    tmp_res <- rbind( tmp1 , tmp2 )
    tmp_res$Gene_Symbol <- gene
    tmp_res$p <- ""
    tmp_res$p <- p
    tmp_res$p_text <- ""
    tmp_res$p_text[1] <- p_text
    result <- rbind( result , tmp_res )
  }
  return(result)
  }
}



plotGene <- function(result=result,y_max=y_max,y_lab=y_lab){
    plot <- ggplot(result, aes(x = Class_use , y = copy_number, color = type , fill = type)) +
    #geom_boxplot(size = 1.2 , outlier.shape = NA ) + ## 去除散点，加粗线
    #geom_jitter(position = position_jitterdodge(0.8) , size = 1) + 
    geom_violin(trim=FALSE) +
    geom_boxplot(width=0.2,position=position_dodge(0.9),fill="white",color="black")+ #绘制箱线图
    #scale_y_log10() +
    facet_wrap(vars(show_text),ncol=7,drop = F) +
    #facet_grid( .~ show_text , scales = "free_x" ) +
    scale_color_manual(values=col) +
    scale_fill_manual(values=col) +
    geom_text(aes(label=p_text , y = y_max , x = 1.5),parse = TRUE,size=5 , color = "black", face='bold') +
    xlab(NULL) +
    ylab(y_lab)+
    theme_bw() +
    theme(
      legend.position = 'none',
      legend.title = element_blank() ,
      panel.grid.major=element_blank(),
      panel.grid.minor=element_blank(),
      panel.background = element_blank(),
      panel.border = element_blank(),
      plot.title = element_text(size = 12,color="black",face='bold'),
      legend.text = element_text(size = 12,color="black",face='bold'),
      axis.text.y = element_text(size = 12,color="black",face='bold'),
      axis.title.x = element_text(size = 12,color="black",face='bold'),
      axis.title.y = element_text(size = 12,color="black",face='bold'),
      axis.text.x = element_text(size = 12,color="black",face='bold') ,
      axis.ticks.length = unit(0.2, "cm") ,
      strip.text.x = element_text(size = 15, colour = "black",face='italic') ,
      axis.line = element_line(size = 0.5)) 
  
  return(plot)
}

###########################################################################################
##导入数据
gistic_IGC = readGistic(gisticAllLesionsFile = gisticAllLesionsFile_IGC,
                        gisticAmpGenesFile = gisticAmpGenesFile_IGC, 
                        gisticDelGenesFile = gisticDelGenesFile_IGC, 
                        gisticScoresFile = gisticScoresFile_IGC, 
                        isTCGA =F)

gistic_DGC = readGistic(gisticAllLesionsFile = gisticAllLesionsFile_DGC,
                        gisticAmpGenesFile = gisticAmpGenesFile_DGC, 
                        gisticDelGenesFile = gisticDelGenesFile_DGC, 
                        gisticScoresFile = gisticScoresFile_DGC, 
                        isTCGA =F)

geneprotein_list <- data.frame(fread(geneprotein_list))
dat_gene1 <- data.frame(fread(all_data_by_genes_IGC))
dat_gene2 <- data.frame(fread(all_data_by_genes_DGC))
colnames(dat_gene1)[1] <- "Gene_Symbol"
colnames(dat_gene2)[1] <- "Gene_Symbol"

dat_info_multi <- data.frame(fread(dat_info_multi))
dat_info_TCGA <- data.frame(fread(dat_info_TCGA))

dat_info_TCGA <- subset(dat_info_TCGA,dat_info_TCGA$From=="TCGA" & dat_info_TCGA$Molecular.subtype=="CIN")
dat_info_TCGA$use_sample <- dat_info_TCGA$Tumor
dat_info_TCGA$use_sample <- gsub("-",".",dat_info_TCGA$use_sample)
dat_info_TCGA$Normal <- dat_info_TCGA$Tumor

dat_info_multi$use_sample <- paste0(dat_info_multi$Tumor,"_",dat_info_multi$Normal)

dat_info <- rbind(dat_info_multi[,c("use_sample","Normal")],dat_info_TCGA[,c("use_sample","Normal")])

###########################################################################################
##数据整合(amp)

gis_all_IGC <- gisticChromGet(gistic = gistic_IGC, fdrCutOff = 0.25, ref.build = "hg19",class="IGC",sub="Amp")
gis_all_DGC <- gisticChromGet(gistic = gistic_DGC, fdrCutOff = 0.25, ref.build = "hg19",class="DGC",sub="Amp")

gis_all_IGC_need <- gis_all_IGC[,c("Cytoband","G_Score","qvalues")]
colnames(gis_all_IGC_need)[2:ncol(gis_all_IGC_need)] <- paste0(colnames(gis_all_IGC_need)[2:ncol(gis_all_IGC_need)],"_IGC")
gis_all_DGC_need <- gis_all_DGC[,c("Cytoband","G_Score","qvalues")]
colnames(gis_all_DGC_need)[2:ncol(gis_all_DGC_need)] <- paste0(colnames(gis_all_DGC_need)[2:ncol(gis_all_DGC_need)],"_DGC")

gis_all <- merge(gis_all_IGC_need,gis_all_DGC_need,by="Cytoband",all = T)

gis_all <- data.frame(sapply(gis_all, function(x) ifelse(is.na(x), 0, x)))

gis_all[, 2:ncol(gis_all)] <- sapply(gis_all[, 2:ncol(gis_all)], function(x) {as.numeric(as.character(x))})

gis_all$Gscore_diff <- gis_all$G_Score_IGC-gis_all$G_Score_DGC
write.table(gis_all,paste0(out_path,"/amp/test.tsv"),sep = "\t",row.names = F,quote = F)
gis_all$logqvalues_IGC <- -log10(gis_all$qvalues_IGC)
dat <- gis_all
q_t <- 0.25
diff_t <- 0.5
sig <- subset(dat, logqvalues_IGC > -log10(q_t) & Gscore_diff > diff_t )
dat$note <- ifelse(dat$logqvalues_IGC > -log10(q_t) & dat$Gscore_diff > diff_t, 'sig', 'non.sig')
dat$logqvalues_IGC <- ifelse(dat$logqvalues_IGC==Inf, 0, dat$logqvalues_IGC)

###########################################################################################
##定位各显著遗传带中的基因
a <- gistic_IGC@data
b <- subset(a,a$Cytoband %in% gis_all_IGC$Unique_Name )
colnames(b)[ncol(b)] <- "Unique_Name"

c <- merge(b,gis_all_IGC,by="Unique_Name")
d <- subset(c,c$Cytoband %in% dat$Cytoband[dat$note=="sig"])

sig$text <- ""
for(i in 1:nrow(sig)){
  sig$text[i] <- paste( sig$Cytoband[i] ,"(" ,length(unique(d$Hugo_Symbol[d$Cytoband==sig$Cytoband[i] ])),")" )
}

write.table(d,paste0(out_path,"/amp/sig_gene_amp_",Class,".tsv"),sep = "\t",row.names = F,quote = F)


###########################################################################################
##IGC/DGC比较

d$Hugo_Symbol <- gsub("^\\[","",d$Hugo_Symbol)
d$Hugo_Symbol <- gsub("\\]$","",d$Hugo_Symbol)
show_symbol <- unique(d$Hugo_Symbol)

result <- datWilcoxTest(dat_gene1=dat_gene1,dat_gene2=dat_gene2,show_symbol=show_symbol)

result$gene_type <- ifelse(result$Gene_Symbol %in% geneprotein_list$gene_name,"coding","noncoding")

outcompare <- result %>%
    group_by( Gene_Symbol , type, p,gene_type ) %>%
    summarize( copy_number = median(copy_number) )
  
out_name <- paste0( out_path , "/amp/CIN_IGCcompareDGC_Driver.copy_number_amp_",Class, ".tsv" )
write.table( outcompare , out_name , row.names = F , quote = F , sep = "\t" )

##基因太多，只显示P<0.05的那些
#result <- subset(result,result$p < 0.05)

###########################################################################################
##画图（类似于火山图）只保留基因显著的遗传带为sig

e <- subset(d,d$Hugo_Symbol %in% unique(result$Gene_Symbol) & d$Hugo_Symbol %in% unique(geneprotein_list$gene_name) )
#sig <- subset(sig, Cytoband %in% unique(e$Cytoband))

sig$text <- ""
for(i in 1:nrow(sig)){
  gene <- unique(e$Hugo_Symbol[e$Cytoband==sig$Cytoband[i]])
  if(length(gene)==0){
    sig$text[i] <- sig$Cytoband[i]
  }else{
    genes <- paste0(gene,collapse=",")
    sig$text[i] <- paste0( sig$Cytoband[i] ,"\n(" ,genes,")")
  }
}

p1 <- plotCytoband(dat=dat,sig=sig,type="Amp")
out_name <- paste0(out_path,"/amp/CIN_IGCcompareDGC_Cytoband_amp_",Class,".pdf")
ggsave(out_name, plot = p1, width = 6, height = 6)

###########################################################################################
##画图（IGC/DGC比较）

col <- c(
  rgb(red=247,green=184,blue=71,alpha=255,max=255) ,
  rgb(red=2,green=100,blue=190,alpha=255,max=255) 
)

names(col) <- c("IGC" , "DGC" )

if(nrow(result) >0){
	sample_class_num <- result %>%
	group_by(type) %>%
	summarize( nums_class = length(unique(Normal)) )

	result <- merge( result , sample_class_num , by = "type")
	result <- subset(result,result$gene_type=="coding")
	result$Class_use <- paste0( result$type , "\n" , "(" , result$nums_class , ")" )
	result$Class_use <- factor( result$Class_use , levels = unique(result$Class_use)[order(unique(result$Class_use) , decreasing=T)] , order = T )

	y_max <- max(result$copy_number) * 1.5
	y_lab <- "log2(copy number/2)"
	result$type <- factor( result$type , levels = c("IGC" , "DGC") )
	## 基于的顺序按照p值由小到达
	result$Gene_Symbol <- factor( result$Gene_Symbol , levels = unique(result[order(result$p),"Gene_Symbol"]))

	## 注释基因位置
	result <- merge( result , dat_gene1[,c("Gene_Symbol" , "Cytoband")] , by = "Gene_Symbol" )
	result$show_text <- paste( result$Gene_Symbol , "\n" , result$Cytoband )
	result$show_text <- gsub( " " , "" , result$show_text)
	result$show_text <- factor( result$show_text , levels = unique(result[order(result$p),"show_text"]))

	sample_class_num <- result %>%
	group_by(type) %>%
	summarize( nums_class = length(unique(Normal)))


	plot <- plotGene(result=result,y_max=y_max,y_lab=y_lab)
	out_name <- paste0( out_path , "/amp/CIN_IGCcompareDGC_Gene_amp_",Class,".pdf" )
	if(length(unique(result$Gene_Symbol))>=7){
	ggsave(file=out_name,plot=plot,width=14,height=length(unique(result$Gene_Symbol))/7 * 3,limitsize = FALSE)
	}else{
	ggsave(file=out_name,plot=plot,width=length(unique(result$Gene_Symbol))/7 * 14,height=3,limitsize = FALSE)
	}
}



########################################################################################################################################
########################################################################################################################################
##数据整合(del)

gis_all_IGC <- gisticChromGet(gistic = gistic_IGC, fdrCutOff = 0.25, ref.build = "hg19",class="IGC",sub="Del")
gis_all_DGC <- gisticChromGet(gistic = gistic_DGC, fdrCutOff = 0.25, ref.build = "hg19",class="DGC",sub="Del")

gis_all_IGC_need <- gis_all_IGC[,c("Cytoband","G_Score","qvalues")]
colnames(gis_all_IGC_need)[2:ncol(gis_all_IGC_need)] <- paste0(colnames(gis_all_IGC_need)[2:ncol(gis_all_IGC_need)],"_IGC")
gis_all_DGC_need <- gis_all_DGC[,c("Cytoband","G_Score","qvalues")]
colnames(gis_all_DGC_need)[2:ncol(gis_all_DGC_need)] <- paste0(colnames(gis_all_DGC_need)[2:ncol(gis_all_DGC_need)],"_DGC")

gis_all <- merge(gis_all_IGC_need,gis_all_DGC_need,by="Cytoband",all = T)

gis_all <- data.frame(sapply(gis_all, function(x) ifelse(is.na(x), 0, x)))

gis_all[, 2:ncol(gis_all)] <- sapply(gis_all[, 2:ncol(gis_all)], function(x) {as.numeric(as.character(x))})

gis_all$Gscore_diff <- gis_all$G_Score_IGC-gis_all$G_Score_DGC
write.table(gis_all,paste0(out_path,"/del/test.tsv"),sep = "\t",row.names = F,quote = F)
gis_all$logqvalues_IGC <- -log10(gis_all$qvalues_IGC)
dat <- gis_all
q_t <- 0.25
diff_t <- 0.5
sig <- subset(dat, logqvalues_IGC > -log10(q_t) & Gscore_diff > diff_t )
dat$note <- ifelse(dat$logqvalues_IGC > -log10(q_t) & dat$Gscore_diff > diff_t, 'sig', 'non.sig')
dat$logqvalues_IGC <- ifelse(dat$logqvalues_IGC==Inf, 0, dat$logqvalues_IGC)

###########################################################################################
##定位各显著遗传带中的基因
a <- gistic_IGC@data
b <- subset(a,a$Cytoband %in% gis_all_IGC$Unique_Name )
colnames(b)[ncol(b)] <- "Unique_Name"

c <- merge(b,gis_all_IGC,by="Unique_Name")
d <- subset(c,c$Cytoband %in% dat$Cytoband[dat$note=="sig"])

write.table(d,paste0(out_path,"/del/sig_gene_del_",Class,".tsv"),sep = "\t",row.names = F,quote = F)


###########################################################################################
##IGC/DGC比较

d$Hugo_Symbol <- gsub("^\\[","",d$Hugo_Symbol)
d$Hugo_Symbol <- gsub("\\]$","",d$Hugo_Symbol)
show_symbol <- unique(d$Hugo_Symbol)

result <- datWilcoxTest(dat_gene1=dat_gene1,dat_gene2=dat_gene2,show_symbol=show_symbol)
result$gene_type <- ifelse(result$Gene_Symbol %in% geneprotein_list$gene_name,"coding","noncoding")

outcompare <- result %>%
    group_by( Gene_Symbol , type, p ) %>%
    summarize( copy_number = median(copy_number) )
  
out_name <- paste0( out_path , "/del/CIN_IGCcompareDGC_Driver.copy_number_del_",Class, ".tsv" )
write.table( outcompare , out_name , row.names = F , quote = F , sep = "\t" )


##基因太多，只显示P<0.05的那些
result <- subset(result,result$p < 0.05)

###########################################################################################
##画图（类似于火山图）只保留基因显著的遗传带为sig

e <- subset(d,d$Hugo_Symbol %in% unique(result$Gene_Symbol) & d$Hugo_Symbol %in% unique(geneprotein_list$gene_name) )
#sig <- subset(sig, Cytoband %in% unique(e$Cytoband))

sig$text <- ""
for(i in 1:nrow(sig)){
  gene <- unique(e$Hugo_Symbol[e$Cytoband==sig$Cytoband[i]])
  if(length(gene)==0){
    sig$text[i] <- sig$Cytoband[i]
  }else{
    genes <- paste0(gene,collapse=",")
    sig$text[i] <- paste0( sig$Cytoband[i] ,"\n(" ,genes,")")
  }
}

p1 <- plotCytoband(dat=dat,sig=sig,type="del")
out_name <- paste0(out_path,"/del/CIN_IGCcompareDGC_Cytoband_del_",Class,".pdf")
ggsave(out_name, plot = p1, width = 6, height = 6)

###########################################################################################
##画图（IGC/DGC比较）

col <- c(
  rgb(red=247,green=184,blue=71,alpha=255,max=255) ,
  rgb(red=2,green=100,blue=190,alpha=255,max=255) 
)

names(col) <- c("IGC" , "DGC" )

if(nrow(result) >0){
	sample_class_num <- result %>%
	group_by(type) %>%
    summarize( nums_class = length(unique(Normal)) )
  
	result <- merge( result , sample_class_num , by = "type")
	result <- subset(result,result$gene_type=="coding")
	result$Class_use <- paste0( result$type , "\n" , "(" , result$nums_class , ")" )
	result$Class_use <- factor( result$Class_use , levels = unique(result$Class_use)[order(unique(result$Class_use) , decreasing=T)] , order = T )

	y_max <- max(result$copy_number) * 1.5
	y_lab <- "log2(copy number/2)"
	result$type <- factor( result$type , levels = c("IGC" , "DGC") )
	## 基于的顺序按照p值由小到达
	result$Gene_Symbol <- factor( result$Gene_Symbol , levels = unique(result[order(result$p),"Gene_Symbol"]))

	## 注释基因位置
	result <- merge( result , dat_gene1[,c("Gene_Symbol" , "Cytoband")] , by = "Gene_Symbol" )
	result$show_text <- paste( result$Gene_Symbol , "\n" , result$Cytoband )
	result$show_text <- gsub( " " , "" , result$show_text)
	result$show_text <- factor( result$show_text , levels = unique(result[order(result$p),"show_text"]))

	sample_class_num <- result %>%
	group_by(type) %>%
	summarize( nums_class = length(unique(Normal)))



	plot <- plotGene(result=result,y_max=y_max,y_lab=y_lab)
	out_name <- paste0( out_path , "/del/CIN_IGCcompareDGC_Gene_del_",Class,".pdf" )
	if(length(unique(result$Gene_Symbol))>=7){
	ggsave(file=out_name,plot=plot,width=14,height=length(unique(result$Gene_Symbol))/7 * 3,limitsize = FALSE)
	}else{
	ggsave(file=out_name,plot=plot,width=length(unique(result$Gene_Symbol))/7 * 14,height=3,limitsize = FALSE)
	}

}

